% jsinusoids
% trying to estimate power spectrum of unevenly spaced data
% to use in psd est of champ data
% latest date 19.11.2003
% uses lombscargle algorithm to estimate power spectrum
% refernce http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=993&objectType=file
function[x,y] = sinu_eef(),

T = 4098*10;
len = 512;
samplig_rate = 5;%minutes


[dummy,f] = pwelch(rand([1,len]),hanning(len),1,len,1/(samplig_rate*60));%just to get frequencies

%used 3 minutes for conveinece 
%will resample at 32 interval (--96 minutes) then the new fft interval if
%512/32=16

fs = 1/(samplig_rate*60);%300seconds of Julia Vz seconds
%x = (0:2879)/fs; %(10 days at 300 secs invl)
%x = (0:19079)/fs;
x = (0:(T-1))/fs; %(10 days at 300 secs invl)

P = log((1./(f*3600)))-4;% -2 is just to reduce amplitude
A = 10.^(P);%to get power increasing 10 times per decade (100 times when using poer spectra)
A(1) = 0;

%y = A(1)*cos(2*pi*x*Frq(1)) + A(2)*sin(2*pi*x*Frq(2)+2*pi/4) + normrnd(0,A(3),1,length(x));
%y = A(1)*sin(2*pi*x*Frq(1)+0)+ A(1)*sin(2*pi*x*Frq(2)+2*pi/4);
% y = A(2)*sin(2*pi*x*Frq(2)+ pi/4) + normrnd(0,A(4),1,length(x));
% note coeff = inversion(x,y,1) just brigs back the power (test for
% inversion) 20 Oct 2005

y = zeros(size(x));
for i = 2:length(f),
    y = y + real(A(i))*cos(2*pi*x*f(i)) + normrnd(0,A(i)/1000,1,length(x));
end;
    
    

